FOSS4G_2022_grass_jupyter slides

Using GRASS GIS in Jupyter Notebooks:
An Introduction to grass.jupyter

Caitlin Haedrich, Vaclav Petras, Anna Petrasova

grass_jupyter

NCSU GeoForAll Lab at the Center for Geospatial Analytics
North Carolina State University

Caitlin Haedrich

PhD student GeoForAll Lab at North Carolina State University
Advised by Helena Mitasova (and Anna Petrasova and Vashek Petras)
Interested in geovisualization and increasing accessibilty of GIS tools
Working with integrating GRASS GIS with Jupyter through Google Summer of Code 2021 and Mini Project Grant from GRASS GIS
TA for graduate-level Geospatial Computation and Simulation

geoforall

ncstate

GRASS GIS

Open-source GIS software with over 500 tools and 400 addons
Interfaces: graphical user interface, command line, Python, C
3rd Party Interfaces: R, QGIS, actinia (REST API)

GUI

grassgis

Project Jupyter

“The Jupyter Notebook is an open-source web application that allows you to create and share documents that contain live code, equations, visualizations and narrative text.” – jupyter.org

Jupyter

In [1]:

import sys

v = sys.version_info
print(f"We are using Python {v.major}.{v.minor}.{v.micro}")

We are using Python 3.9.5

Teaching Geospatial Analytics with Jupyter

GIS714: geospatial computation and simulation: graduate level course required for Geospatial Analytics PhD students

oldGIS714

grass.jupyter

Released with main GRASS GIS distribution starting in version 8.2
Python package that makes GRASS more accessible in Jupyter by providing:
- Data management functions
- Visualizations classes for viewing and interacting with GRASS through in-line figures
Created with syntax consistent with GRASS's existing Python API and command line interface

You can also try it out in Binder! (also linked from the GRASS GIS GitHub README)

Getting Started

In [4]:

# Import Python standard library and IPython packages we need.
import subprocess
import sys

# Ask GRASS GIS where its Python packages are.
sys.path.append(
    subprocess.check_output(["grass80", "--config", "python_path"], text=True, shell=True).strip()
)

# Import the GRASS GIS packages we need.
import grass.script as gs
import grass.jupyter as gj

# Start GRASS Session
session = gj.init("data", "nc_spm_08_grass7", "PERMANENT")

# Make a new mapset for this assignment
gs.run_command("g.mapset", mapset="HW3_water_simulation", location="nc_spm_08_grass7", flags="c")

# Set region to the high resolution study area
gs.run_command("g.region", region="rural_1m")

And then we are ready to begin using GRASS!

Since grass.jupyter is written in Python, we'll use the GRASS Python API.

Teaching Geospatial Computation and Simulation with grass.jupyter

GIS714: Geospatial Computation and Simulation, spring 2022,

Course Repo on GitHub

We'll use an example from the course on water simulation...

Create a map of flooded area¶

We create a map of flooded area with r.lake (GRASS manual for r.lake) by providing a water level and a seed point:

In [21]:

gs.run_command("r.lake", elevation="elev_lid792_1m", water_level=113.5, lake="flood1", coordinates="638728,220278")

# See results
flood1 = gj.Map(use_region=True)
flood1.d_rast(map="elev_lid792_1m")
flood1.d_rast(map="flood1")
# Display map
flood1.show()

Out[21]:

Question 4¶

Increase water level to 113.7m and 114.0m and create flooded area maps at these two levels.

In [22]:

#### Your Answer Here

In [ ]:

# 3D map example

grass.jupyter Syntax

Calling display modules using Python magic short-cut:

d.rast -> Map.d_rast(), Map3D.d_rast

Estimating inundation extent using HAND methodology¶

We will use two GRASS addons, r.stream.distance and r.lake.series, to estimate inundation with Height Above Nearest Drainage methodology (A.D. Nobre, 2011). We need to install them first:

In [ ]:

#!g.extension r.stream.distance
#!g.extension r.lake.series

For this section, we will change our computation region to elevation which is a larger study area than we used above. Because we are using a new region (and also a higher threshold of 100,000), we need to run r.watershed again to compute the flow accumulation, drainage and streams. We convert the streams to vector for better visualization.

In [24]:

gs.run_command("g.region", raster="elevation")
gs.run_command("r.watershed", elevation="elevation", accumulation="flowacc", drainage="drainage", stream="streams_100k", threshold=100000)
gs.run_command("r.to.vect", input="streams_100k", output="streams_100k", type="line")

Now we use r.stream.distance with output parameter difference to compute new raster where each cell is the elevation difference between the cell and the the cell on the stream where the cell drains. This is our HAND terrain model.

In [ ]:

gs.run_command("r.stream.distance", stream_rast="streams_100k", direction="drainage", elevation="elevation", method="downstream", difference="above_stream")

With r.lake.series, we can create a series of inundation maps with rising water levels. r.lake.series creates a space-time dataset. We can use temporal modules to further work with the data...

In [ ]:

gs.run_command("r.lake.series", elevation="above_stream", start_water_level=0, end_water_level=5, water_level_step=0.5, 
               output="inundation", seed_raster="streams")

... or visualize the flood with TimeSeriesMap.

In [ ]:

flood_map = gj.TimeSeriesMap(use_region=True)
flood_map.add_raster_series("inundation")
flood_map.d_legend(color="black", at=(10,40,2,6)) #Add legend
flood_map.d_barscale()
flood.show()

folium Integration

folium...

is a popular Python library for creating Leaflet.js maps

creates interactive HTML maps that can be displayed in-line

has built-in tilesets from OpenStreetMap, Mapbox, and Stamen, and supports custom tilesets

foliumLogo

Moving data from GRASS GIS location projection to folium is a challenge

folium is projected in Web Mercator (EPSG:3857)

However, any coordinates (i.e. vector data) need to be specified in degrees of latitude and longitude (WGS84, EPSG:4326)

folium reprojects latitude and longitude to Web Mercator internally

We pass data to folium by reprojecting to temporary locations

grass_folium

In [ ]:

import folium 

# Create figure
fig = folium.Figure(width=600, height=400)

# Create a map to add to the figure later
m = folium.Map(tiles="Stamen Terrain", location=[35.761168,-78.668271], zoom_start=13)

# Create and add elevation layer to map
gj.Raster("elevation", opacity=0.5).add_to(m)

# Do some cool folium stuff!
# Like make a tooltip
tooltip = "Click me!"
# and add a marker
folium.Marker(
    [35.781608,-78.675800], popup="<i>Center For Geospatial Analytics</i>", tooltip=tooltip
).add_to(m)

# and a circle
folium.Circle(
    radius=120,
    location=[35.769781,-78.663160],
    popup="Great Picnic Area",
    color="crimson",
    fill=False,
    tooltip=tooltip
).add_to(m)

# Add the map to the figure
fig.add_child(m)

In [7]:

# Display figure
fig

Out[7]:

The InteractiveMap class allows users unfamiliar with folium to produce maps easily.

In [8]:

# Create Interactive Map
fig = gj.InteractiveMap(width = 600)
# Add raster, vector and layer control to map
fig.add_raster("elevation", title="Elevation Raster", opacity=0.8)
#fig.add_vector("roadsmajor")
fig.add_layer_control(position = "bottomright")
# Display map
fig.show()

Out[8]:

Course Reflections

Notebook format was generally well received:

Several students noted they enjoyed being able to experiment with GRASS tools by varying parameters and visually observing the effect while working through the tutorial sections of the notebooks.
By hosting the notebooks on GitHub and including a link to Binder, students could easily access and run the notebooks without installing and configuring software on their personal computers.
- This also gave the students exposure to working with Git and GitHub, an important tool for open science and collaborative research.

There were several areas of the notebook deployment that needed improvement or that posed a difficulty to the class:

First, running the notebooks locally was difficult for many students. Due to the required configuration of path variables on Windows operating system and installation of required packages, many students switched to Binder right away.
For students that did get the notebooks running locally, we often encountered issues that were specific to their version of GRASS GIS or operating system.
We also found that Binder was not a reliable host; occasionally, there were too many users accessing Binder’s computational nodes and the notebooks failed. In the future, we may use a JupyterHub hosted by NC State University to circumvent such technical issues.
Secondly, the notebooks assume users have basic Python abilities, are familiar with the fundamentals of geospatial analysis (such as common data formats and projections) and have some experience using Jupyter Notebooks. For students without any of these prerequisite skills, the course could be challenging even with the removal of the command line interface and GUI workflows present in the previous tutorial format.